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o : 

Superconductivity in intercalated graphite CaC6 and H under extreme pressure, in the framework 
of superconducting density functional theory, is discussed. A detailed analysis on how the electron- 
phonon and electron-electron interactions combine together to determine the superconducting gap 
and critical temperature of these systems is presented. In particular, we discuss the effect on the 
calculated T c of the anisotropy of the electron-phonon interaction and of the different approximations 
^f*j ' for screening the Coulomb repulsion. These results contribute to the understanding of multigap and 

anisotropic superconductivity, which has received a lot of attention since the discovery of MgE>2, 
and show how it is possible to describe the superconducting properties of real materials on a fully 
ab-initio basis. 
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V ■ I. INTRODUCTION 

u 

The discovery of two-gap superconductivity in MgB2 has triggered a huge revival of interest on the search of new 
electron-phonon superconductors. The recent experience in the vast majority of electron-phonon superconductors 
has shown that ab — initio calculations can provide a detailed and useful description of the normal state of these 
systems, including dynamical properties and electron-phonon coupling. The treatment of the superconducting state, 
however, is a much more difficult task. On the theoretical side, most of nowadays calculations are based on the 
Eliashberg theory [3, 0,0). While in the latter the electron-phonon interaction is perfectly accounted for, the effects of 
the electron-electron Coulomb repulsion are condensed in a single parameter, fj,* , which is difficult to calculate from 
first principles and which, in most practical applications, is treated as an adjustable parameter usually fitted to the 
O . experimental T c . In this sense the Eliashberg theory, in spite of its tremendous success, has to be considered as a 
semi-phenomenological theory. However, the possibility to describe the superconducting state on a fully ab — initio 
ground and, as a consequence, to be able to predict the superconducting properties of materials, is still highly sought 
> for. 

The density functional theory (DFT), a very successful standard approach in normal state electronic structure 
calculations, has been recently extended to deal with the superconducting state (SCDFT) [27t [28j. SCDFT contains 
no adjustable parameters, and the final critical temperature is the result of material specific quantities, all of them 
computed ab — initio. SCDFT is able to treat superconductors with a wide range of couplings, as shown by several 
investigations, MgB2 being one of the most non-trivial one [4j. 

Superconductivity is the result of a subtle competition between two opposite effects: the phonon mediated attraction 
("e-ph" in the following) and the direct Coulomb repulsion ("e-e") between the electrons. In this work, we show 
how different approximations for the e-e and e-ph interactions (e.g. an isotropic approximation, or an insufficiently 
elaborated screening entering the e-e term) will depart from the experimental results. Considering intercalated 
graphite CaCg and H under pressure as test cases, we will quantify, without ad-hoc assumptions, the effect on T c of 
the superconducting gap anisotropy and of the assumed screened Coulomb potential. The choice of these two systems 
comes from theirpeculiarities: CaCg is a superconductor with a moderate degree of anisotropy [hJ, recently confirmed 
by experiments 0, @|, with electronic bands having different orbital characters. With a critical temperature as high as 
11.5 K, CaC6 exhibits the highest T c [7| among graphite intercalated compounds. DFT calculations of CaC6 electronic 
and dynamical properties pointed out that the e-ph coupling is large enough to yield the observed T c . In CaC6, 
the phonons mostly contributing to superconductivity belong to an optical branch involving Ca displacements and, to 
a lower extent, to two C- related branches at much higher frequency. The different phonon branches couple to different 
Fermi surface (FS) sheets, as pointed out in Ref. QiJ 

High temperature superconductivity in H was suggested 40 years ago by Ashcroft[l2| and has been the subject of 
several investigations (see Ref. [ill and references therein). Unlike CaCg, molecular hydrogen under pressure involves 
basically only one type of atomic orbital, it has very large phonon frequencies (due to the light H mass) and large 
e-ph coupling fill ] . 

The paper is organized as follows: In Section [TT] we summarize the main features of SCDFT and describe our 
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computational approach; in section IIIII we discuss the different approximations used for the Coulomb interaction; in 
Sections IIVI and [V] we present our results for CaCg and H, respectively; finally, in Section IVIi we summarize our 
conclusions. 



II. THE DENSITY FUNCTIONAL THEORY FOR SUPERCONDUCTORS 



The density functional theory [171 ] has enjoyed increasing popularity as a reliable and relatively inexpensive tool to 
describe real materials. In this section we will briefly outline the DFT approach to superconductivity, and refer to the 
original papers for more details. In order to give an introduction to SCDFT, it is instructive to recall how magnetism 
is treated within the DFT. The Hohenberg-Kohn (HK) theorem [l8[ states that all observables, in particular also 
the magnetization, are functionals of the electronic density alone. This, however, assumes the knowledge of the 
magnetization as a functional of the density. Finding an approximation for this functional is extremely hard and, in 
practice, one chooses a different approach. The task can be vastly simplified by treating the magnetization density 
m(r), i.e., the order parameter of the magnetic state, as an additional fundamental density in the density functional 
framework [l9| . An auxiliary field - here a magnetic field B oxt (r) - is introduced, which couples to m(r) and breaks 
the corresponding (rotational) symmetry of the Hamiltonian. This field drives the system into the ordered state. If 
the system is actually magnetic, the order parameter will survive when the auxiliary perturbation is quenched. In 
this way, the ground-state magnetization density is determined by minimizing the total energy functional (free energy 
functional for finite temperature calculations) with respect to both the normal density and the magnetization density. 
Within this approach a much simpler approximations to the xc functional (now a functional of two densities) can lead 
to satisfactory results. 

The same idea is also at the heart of density functional theory for superconductors, as formulated by Oliveira, Gross 
and Kohn [l6[. Here the order parameter is the so-called anomalous density, 

X (r,r') = (* T (r)* i (r')), (1) 

and the corresponding potential is the non-local pairing potential A(r, r'). It can be interpreted as an external pairing 
field, induced by an adjacent superconductor via the proximity effect. Again, this external field only acts to break 
the symmetry (here the gauge symmetry) of the system, and is quenched at the end of the calculation. As in the case 
of magnetism, if the system is actually a superconductor the order parameter will be sustained by the self-consistent 
effective pairing field. The approach outlined so far captures, in principle, all the electronic degrees of freedom. To 
describe conventional phonon-mediated superconductors, also the electron-phonon interaction has to be taken into 
account. 

In order to treat both weak and strong electron-phonon coupling, the electronic and the nuclear degrees of freedom 
have to be treated on equal foo ting . This can be achieved by a multi-component DFT, based on both the electronic 
density and the nuclear density [2Ct l2l|. In addition to the normal and anomalous electronic densities, we also include 
the diagonal of the nuclear density matrix [44J 



r(R) = ($t(Ri) . . . &{R N )$(R N ) . . . $(R X )), (2) 

where $(R) is a nuclear field operator. 

In order to formulate a Hohenberg-Kohn theorem for this system, we introduce a set of three potentials, which 
couple to the three densities described above. Since the electron-nuclear interaction, which in conventional DFT 
constitutes the external potential, is treated explicitly in this formalism, it is not part of the external potential. The 
nuclear Coulomb interaction U nn already has the form of an external many-body potential, coupling to T(R), and 
for the sake of the Hohenberg-Kohn theorem, this potential will be allowed to take the form of an arbitrary N-body 
potential. All three external potentials are merely mathematical devices, required to formulate a Hohenberg-Kohn 
theorem. At the end of the derivation, the external electronic and pairing potentials will be set to zero while the 
external nuclear many-body potential to the nuclear Coulomb interaction. 

As usual, the Hohenberg-Kohn theorem guarantees a one-to-one mapping between the set of the densities 
{n(r), %(r, r'), r(R)} in thermal equilibrium and the set of their conjugate potentials {i£ xt (r) — /z, A ext (r, r'), t£ xt (R)}. 
Therefore all the observables are functionals of the set of densities. Finally, it assures that the grand canonical po- 
tential, 

il[n,x,r]=F[n, X ,r} + |d 3 rn(r)[< xt (r)- M ]-|d 3 r|dV [ X (r, r')A c * xt (r, r ') + h.c] + Jd 3 R r(R)< xt (R), (3) 
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is minimized by the equilibrium densities. We use the notation A[f] to denote that A is a functional of /. The 
functional F[n,\,T] is universal, in the sense that it does not depend on the external potentials, and is defined by 

F[n, x , r] = r>, x , r] + r>, x , r] + t/ c >, x, r] + c/ c >, x , r] - ±s[n, x, r] , (4) 

where S is the entropy of the system, 

S[n, x, T] = -Tr{/3 [n, x, T] ln(po[n, x, T])} . (5) 

In standard DFT one normally defines a Kohn-Sham system, i.e., a non-interacting system chosen such that it has 
the same ground-state density as the interacting one. The variational procedure for this system gives Schrodinger-likc 
(Kohn-Sham) equations for non-interacting electrons subject to an effective (Kohn-Sham) potential. These equations 
are nowadays routinely solved by solid state theorists. In our formalism, the Kohn-Sham system consists of non- 
interacting (superconducting) electrons, and interacting nuclei. We will not describe here the details of the method, 
and will only outline its basic features: The Kohn-Sham potentials, which are derived in analogy to normal DFT, 
include the external fields, Hartree, and exchange-correlation terms. The latter account for all many-body effects of the 
electron-electron and electron-nuclear interactions. Obtaining their explicit form has represented a major theoretical 
effort[22|, [H, [24[. Once this problem has been solved, the problem of minimizing the Kohn-Sham grand canonical 
potential can be transformed into a set of three differential equations that have to be solved self-consistently: One 
equation for the nuclei, which resembles the familiar nuclear Born-Oppenheimer equation, and two coupled equations 
which describe the electronic degrees of freedom and have the algebraic structure of the Bogoliubov-de Gennes [25j 
equations. 

The resulting Kohn-Sham Bogoliubov-de Gennes (KS-BdG) equations read (we use atomic Rydberg units) 



' 2 

Yl 
2 



+ t£(r)-/i 



Mnk(r) + Jd 3 r' A s (r, r')v„ k (r') = E nk u nk {r) , (6a) 
Unk(r) + J dV A s *(r,r> rlk (r') = £„ k w„ k (r) , (6b) 



where w„k(r) and «nk( r ) are the particle and hole amplitudes. This equation is very similar to the Kohn-Sham 
equations in the OGK formalism [lq . However, in the present formulation the lattice potential is not considered 
an external potential but enters via the electron-ion Hartree term. Furthermore, our exchange-correlation potentials 
depend on the nuclear density matrix, and therefore on the phonons. Although Eq.([(j]) and the corresponding equation 
for the nuclei have the structure of static mean-field equations, they contain, in principle, all correlation and retardation 
effects through the exchange-correlation potentials. 

These KS-BdG equations can be simplified by the so-called decoupling approximation [2(1 H3| , which corresponds 
to the following ansatz for the particle and hole amplitudes: 

Unk(r) ~ M„kVnk(r) ; "nk(r) ~ Wnk^nkM , (7) 

where the wave functions y n k(r) are the solutions of the normal Schrodinger equation. In this way the eigenvalues in 
Eqs. © become E n ^ = ±S„k, where 



£„k = V^k + |A nk | 2 , (8) 

and £„k = e n k~ /■*• This form of the eigenenergies allows us to interpret the pair potential A„k as the gap function of the 
superconductor. Furthermore, the coefficients u„k and i^k are given by simple expressions within this approximation 

u nk = -i=sgn(£; n k)e i0 " k Jl + £± , (9a) 



V2 




Finally, the matrix elements A„k are defined as 

A„ k = fd 3 r f dV ^,* ik (r)A s (r, r> nk (r') , (10) 
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and <j> n k is the phase e 1 ^™ 11 = A n k/|A n k|- The normal and the anomalous densities can then be easily obtained from: 



l-l^tanhf^nk 



knk(r)| 2 (11a) 



X (r, r') = \Y, tanh f f ^ ) ^ k (r)^ k (r') . (lib) 

nk 



Within the decoupling approximation, we finally arrive at an equation for the k-resolved superconducting gap A„k, 
which has the following formfH, H3, HI, 



tanh 



1 'ami 15 &„'k' y 
A„k = — -Z„kA„k — - ^nk.n'k' ^ A„<k' • (12) 



^ „ , -C'n'k 
n k 

Eq. ([12"]) is the central equation of the DFT for superconductors. The kernel JC consists of two contributions 
JC = /C c ~ ph + IC°~ C , representing the effects of the e-ph and of the e-e interactions, respectively. The diagonal term 
Z plays a similar role as the renormalization term in the Eliashberg equations. Explicit expressions of /C c_ph and Z, 
which are the results of the approximate functionals, are given in Eqs. 9 and 11 of Ref. [28| respectively. These two 
terms involve the e-ph coupling matrix, while JC°~ C contains the matrix elements of the screened Coulomb interaction 
(the explicit expression is given below). Eq. lfl2|) has the same structure as the BCS gap equation, with the kernel 
JC replacing the model interaction of BCS theory. This similarity allows us to interpret the kernel as an effective 
interaction responsible for the binding of the Cooper pairs. Moreover, we emphasize that Eq. (|12p is not a mean-field 
equation (as in BCS theory), since it contains correlation effects via the SC exchange-correlation functional entering 
JC and Z. Furthermore, it has the form of a static equation - i.e., it does not depend explicitly on the frequency - and 
therefore has a simpler structure (and computationally more manageable) than the Eliashberg equations. However, 
this certainly does not imply that retardation effects are absent from the theory. Once again, retardation effects enter 
through the xc functional, as explained in Refs. [U [28| . The SCDFT allows to treat on the same footing the e-ph 
and the screened e-e interactions. These terms, however, can be treated at different levels of approximation. 

We calculated the screened Coulomb matrix elements (ME) with respect to the Bloch functions, for the whole 
energy range of relevant valence and conduction states. The different nature of the electronic bands in each material 
(e.g. some of them can be highly localized while others more delocalized), strongly calls for the use of a non-diagonal 
screening, including local field effects. In order to properly describe these effects, a very important step to achieve 
good agreement with the experiment, we calculated the static random phase approximation (RPA) dielectric matrix 
(DM) e^ 1 (q,G,G') 1 using the pseudopotential-based SELF code (3(j. The explicit expression of the kernel JC e ~ e in 
reciprocal space reads: 

(n / fcyC«+ G > ,r |nJfe) lnk\e- l( y q+G ')- r \n'k') 

SFTWgi (13) 

where q = k' — k. The detailed structure, energy and k-dependence of the kernel JC e ~ e changes the Coulomb 
renormalization effect (due to the different scales of the vibrational and electronic energies) considered: In many 
materials, it only acts as a scaling of the superconducting gap at the Fermi energy, thus keeping the main structures 
of the coupling (See Ref. [ToL [3ll.l32l.l33j). On the other hand, in strongly anisotropic materials with small interband 
interaction (like MgB2), the Coulomb renormalization turns out to be non trivial 0, [34|. 

In this paper we report a study of the role of the different approximations for the Coulomb and the electron-phonon 
interactions on the solution of the SCDFT gap equation. This is illustrated in CaC6[10( and H under extreme pressure. 
As mentioned above, the normal state calculations, necessary for the study of the superconducting state, are performed 
within DFT in the LDA or GGA approximations. Computationally, the electronic and dynamical properties are 
obtained using the pseudopotential method as implemented in the QUANTUM-ESPRESSO package [35j: the screened 
Coulomb matrix elements are obtained with the SELFI30| code. 



III. EVALUATION OF COULOMB MATRIX ELEMENTS 



In this section we discuss the different approximate formulations of the Coulomb interaction analyzed in this work. 
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A. Sham-Kohn 



The first and simplest approximation includes a Thomas- Fermi (TF) dielectric function together with free-electron 
wavef unctions. This leads to isotropic (i.e. k independent) matrix elements, expressed by the analytic formula [27l. |28| : 



Ktlf = — = 



In 



£ + C + 2y^r+^ F /2 



(14) 



where the Fermi level is fixed by the number of valence electrons and the TF wavevector is that of a jellium having the 
same density of states at the Fermi energy of the true system N(0), i.e. q^ F = 87riV(0) . A justification for the use of 
this approximation, in the spirit of the seminal paper of Sham and Kohn |39j|. is given in Ref. [27| . Typically, Eq. (|14|) 
provides reasonable results for systems with delocalized electrons [28l [33 ] . Furthermore, it avoids the cumbersome 
calculation of the anisotropic Coulomb matrix elements, calculated with the Kohn-Sham Bloch functions of the system 
(see below). 



B. Bloch wavefunctions and Thomas-Fermi 



In order to include the effects of the wavefunction localization in the evaluation of the matrix elements, we used 
the real Kohn-Sham Bloch states of Eq. (fl~3|) . together with the TF screened dielectric function. The TF vector is 
again q^ F — 8irN(0), where N(0) is the density of states at the Fermi energy in the real material. In this way, we 
correct for the locality of the wavefunctions in the evaluation of the Coulomb integrals, still keeping, for the screening, 
a very approximate form. The present approximation and the ones described below in this section correspond to an 
anisotropic treatment of the Coulomb interaction. 



C. Diagonal random phase approximation (RPA) 



The TF dielectric function is only valid in the limit q — > 0. To improve this aspect we calculate the dielectric 
function e(q, G, G") in the RPA, initially limiting ourselves to the diagonal form, where we set to zero all the G^G' 
terms. In real space this corresponds to a screening depending not on r and r' separately but only onr-r'. On the 
other hand, the electronic polarizability is built by taking into account the real band structure of the system. This 
approximation is expected to be reasonable for closely packed systems with delocalized electrons. 



D. Non-diagonal RPA 



With a fully non local RPA screening (with non-zero G ^ G' terms) the electronic response will depend on the two 
spatial indices r and r' separately. This improves the efficiency of the screening for strongly localized states, with a 
resulting reduction of the Coulomb matrix elements in the corresponding regions. As an example, this improvement 
affects mainly the graphene sheets in CaCg. 



IV. RESULTS ON CaC 6 



Before discussing the various approximations for the Coulomb interaction, we briefly describe the electronic states 
of CaCg. This material crystallizes in the rhombohedral space group R3m [7j, and has one unit formula per primitive 
cell. The rhombohedral lattice parameter is 5.17 A with rhombohedral angle 49.55°. We plot in Fig. Q] the band 
structure of the system. Three bands cross the Fermi energy. Two of them have p character and come from the tt 
bonding pattern of the graphene layer. They give rise to two 2-dimensional Fermi surfaces (Fig.[TJ right panel). The 
third band corresponds to a spherical Fermi surface with the charge localized in the interlayer space and on the Ca 
ion. Its character is mainly s (d) in the k x , k v (k z ) directions (z-axis being normal to the graphene layer). 

The presence of an intercalant band crossing the Fermi energy has been shown to be a necessary condition for 
superconductivity in graphite intercalated compounds (GICs), as pointed out by Csanyi et al. (36|. Interlayer bands, 
uncorrelated to a particular atomic orbital nature, have been largely discussed in graphite GICs [§■ 1361 l37l l3q and 
MgB 2 i.4]. In CaCg, however, this band gains a further strong contribution from the Ca 4s and 3d orbitals. This band 
is coupled to the in-plane Ca optical phonon branch; these low frequency modes provide a strong contribution to the 
e-ph constant A 0, 0, fic| . The energy window between -20 and -4 eV is mainly occupied by the bands forming the 
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l r % x r to 5 

DOS 



FIG. 1: Left panel: Electronic band structure of CaC6, with the Fermi energy set to zero. Different 
color for the Kohn-Sham eigenvalues are used corresponding to the s (black) p (orange-light grey) or 
d (red-gray) main orbital character. Central panel: Total density of states (DOS) (thick black line), 
projected DOS on Ca site (dashed red-gray line) and on C states, divided in sp 2 part (turquoise-light 
gray), and p z part (long dashed black). Right panel: CaCe Fermi surface, showing the external 7r 
sheet (1, green), the Ca-related sphere (2, blue) and the internal tt sheet (3, red). 



sp 2 bonds pattern of graphene, as shown by the orbital decomposed density of states in Fig. [TJ The p z bands starts 
around the L symmetry point (corresponding to the A point in the hexagonal unit cell, i.e. to the zone boundary in z 
direction) at about -10 eV. Just above the Fermi level, mainly between 2 and 5 eV the DOS shows a peak originated 
by a set of flat bands due to Ca d states. Above 5 eV the percentage of the total charge projected over atomic orbitals 
starts to reduce corresponding to the free electron limit of the KS-wavefunctions. 

A. Anisotropic approach for the Coulomb interaction 

We now discuss the effect of the approximate treatments of the Coulomb interaction described in Sec. IIIII on the 
superconducting T c . In the following analysis, we assume a fully anisotropic phononic kernel. Within the Sham-Kohn 
approximation we obtain T c =13.5 K, to be compared with the experimental value T° xp =11.5 K. The Sham-Kohn 
approximation is rather simple and gives quite reliable results for nearly free electron metals [28l. l32l|. However, for 
CaC6 this approximation is not expected to hold, since it neglects the presence of localized states with much stronger 
Coulomb interactions. Therefore in the Sham-Kohn approximation an overestimate of the superconducting gap and 
T c is expected and found. 

With the Thomas-Fermi approximation and Bloch wavefunctions we obtain T c rs 8.2 K, underestimating the 
experimental value. Remarkably, the critical temperature obtained using the RPA approximation with diagonal 
screening turns out to be T c « 8.0 K, very close to the one found in the TF approach and quite smaller than the 
experimental value. This result shows that the TF screening is quite similar to the static diagonal RPA one, and that 
indeed the T c underestimation is related to the lack of local field corrections. A diagonal screening in fact averages 
out over a unit cell of the system. A correct screening, on the other hand, should be more effective in regions of the 
unit cell with larger charge density. As a consequence, inclusion of non-diagonal terms in the screening produces an 
increased T c of 9.4 K. This value, while not in excellent agreement with the experiment, represents an improvement 
over the previous T c values (we recall that in our calculations no fitting parameter to the experiment is used) . In the 
Table below we resume the superconducting properties in the discussed cases, including T c and the value of the gap 
A„fc on the different Fermi surface sheets. 
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SK -TF 


wf -TF 


diag - RPA 


RPA 


no Coulomb 


(Ai) [meV] 


2.20 


1.24 


1.19 


1.42 


5.40 


(A 2 ) [meV] 


2.94 


1.78 


1.77 


2.00 


6.72 


(As) [meV] 


2.49 


1.39 


1.37 


1.67 


6.08 


T C [A1 


13.5 


8.2 


8.0 


9.4 


31 



The brackets < • ■ • > indicate a Fermi surface average and the indices 1,2,3 correspond to the three parts of the Fermi 
surface (see Fig. [TJ, namely the external (green) ir surface, the Ca sphere (blue) and the internal (red) ir surface. 

The detailed behavior of the gap A„k as a function of the normal state eigenvalue is given in Fig. O The set of 
black points shows the gap at T = K, as a function of the energy distance from the Fermi level Ep. For each energy 
the gap is not a single-valued function, i.e. the gap is anisotropic in the reciprocal space. In particular, the gap values 
exhibit a moderate anisotropy, rather than a multigap character. 
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FIG. 2: CaC6 superconducting gap A n k as a function of the energy distance from the 
Fermi energy (n in the figure). The right panel shows the energy distribution (in arbitrary 
units) of the gap, resolved over the three Fermi surfaces. The histogram is evaluated 
within the energy window from -2 to 2 meV. 



B. Isotropic approach for the Coulomb interaction 



All the above results (excluding the Sham-Kohn approximation) refer to fully anisotropic calculations, in which 
the k, k' dependence of the Coulomb matrix elements is included in the solution of the SCDFT gap equation. We 
also calculated T c with an isotropic approach in which the K^ikn'k' matrix elements are averaged in k and k' over 
isoenergy surfaces, yielding a corresponding V (£, £') function. The average is done for each one of the approximations 
discussed. The critical temperatures turn out to be quite close (of the order of few percent) to the corresponding fully 
anisotropic ones, showing that the anisotropy effect on the Coulomb repulsion is almost negligible, unlike in the case 
of MgB2 [1, H3| . On the other hand (see Sec. IIV C[) . this is not true for the phononic part of the kernel [al [til. [Io| . 

A graphical comparison between the (averaged) different approximations is given in Fig. [31 where we report (upper 
panel) two cuts of the function V namely the diagonal part V (£,£) and the elements with the first energy 

corresponding to the Fermi level V (0, £'); in the lower panel a graphical representation of the whole function in the 
TF and in RPA approximation is presented. The Sham-Kohn approach gives a smooth function, with no correlation 
with the electronic structure of the material. The upturn at about -18 eV corresponds to the energy of the k = 
state in a free electron gas with the same number of electrons per unit cell. The use of TF approximation with the 
Bloch wavefunctions shows the effect of the orbital nature of states reflected in the matrix elements. In Fig. [3] for 
example we see the strong Coulomb repulsion between very localized Ca-d states giving a peak in the Coulomb matrix 
elements for £ ~ £' ~ 4 eV. The strong localization implies a nearly zero interaction between these states and the rest 
of the electronic structure (the darkest regions in the matrix), whose charge is located in different areas . Moreover 
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-20 -10 10 1 2 3 4 5 6 

Energy [eV] DOS [states /(eV cell)] 

FIG. 3: Comparison between different approximations in calculating the Coulomb matrix elements 
(ME). The central panel shows, with a scale of colors, the matrix of averaged Coulomb ME (defined 
in Eq. 1 13 p . Because this matrix is symmetric, we used the upper-left part to show Thomas- Fermi 
ME (using Bloch wavefunctions) and the bottom-right part for RPA ME. The upper panel shows 
two cuts of this matrix for the various approximations discussed in the text: A diagonal cut [V (£,£), 
straight lines] and an horizontal cut passing for Ef [V (0, £'), dashed lines]. On the right the density 
of electronic states. 



we can distinguish the signature of the sp 2 bonds in the region -20 to -5 eV in a sharp increase of the matrix elements 
with poor interaction with higher energy states. 

Going from Thomas-Fermi to diagonal-RPA screening does not introduce any important changes. Local fields 
effects, on the other hand, strongly reduce the matrix elements below the Fermi level. In the high energy region 
all the approximations behave in a similar way as the states are free-electron like. The RPA including local field 
effects, although giving similar structures with respect to the Bloch wavefunctions and TF approach, introduces a 
more efficient screening in the graphene layer so that the corresponding matrix elements are reduced of about 20%. 

C. Effects of anisotropy in the phonon coupling 

It is interesting, at this point of the discussion, to analyze the role of the e-ph interaction anisotropy and its effects 
on T c and on the anisotropic character of the superconducting gap. 

From a multiband BCS model (i.e. the Suhl-Matthias- Walker model [Ilj]), it is known that an anisotropy in the 
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phonon coupling always yields a T c enhancement with respect to an averaged coupling. To discuss this kind of 
anisotropy, it is useful to compare CaCg with MgB 2 , the prototypical two gap superconductor. We first define the 
band resolved e-ph coupling Ay = / a 2 l Fij(uj) / 'ivdui and Eliashberg functions a 2 Fij 



(15) 



k,k' v 



where: i and j are indices of the selected sheet of the Fermi surface; Nj is partial DOS at Fermi level, for the sheet 
are the e-ph matrix elements; v is the phonon branch; cj q „ the phonon frequencies. Moreover, we define 



the DOS-rcnormalizcd, band resolved matrix elements V? A = A 



In the following, we present the partial DOS and couplings for CaCg and MgB2: For CaC6 we separate the Fermi 
surface in three parts: (1) external n bands, (2) Ca sphere, (3) internal tt bands (the one cutting the Ca sphere, as 
shown in Fig.[T|). We obtained: 



A : 

0.286 0.173 0.223 
0.518 0.315 0.425 
0.382 0.245 0.303 



V ph [eV] : 

0.363 0.651 0.492 
0.651 1.177 0.932 
0.492 0.932 0.678 



DOS 



eV-ccllJ 



0.79 
0.27 
0.45 



For MgB 2 we have [34 



A : 

1.017 0.213 
0.156 0.448 



VP' 1 [eV] : 

3.390 0.520 
0.520 1.093 



DOS [ 



eV-ccllJ 



0.30 
0.41 



where the first entry in each row or column corresponds to the two a bands and the second to the it [42j . 

These data show clearly that both systems are quite anisotropic: In MgB 2 , the much stronger coupling in the a 
rather than in the tt Fermi surfaces and the small interband coupling, makes the calculated anisotropic T c double 
with respect to the isotropic one 0, [Hj]. Within the Suhl-Matthias- Walker model [4l[ this can be roughly explained 
by a large X m ax = 1-070 (defined as the maximum eigenvalue of the Ay matrix - very similar to An) with respect to 
the average, isotropic coupling A=0.860. 

In CaC6, instead, both Ca and C bands contribute strongly to the global coupling. Consider the Ca spherical FS, 
that in CaC6 has the largest superconducting gap. As shown by the Ay matrix, the contribution coming from the 
interband scattering is larger than the intraband one (A2i=0.518 and A 2 3 = 0.42 5, against the intraband A 2 2=0.315): 
Although the Ca intraband average coupling (as represented by the V? matrix) is the largest one, the large phase 
space (DOS) for interband scattering, with a reasonably large interband scattering potential, makes the difference. 
We calculated A ma :r=0.912 and A=0.855, both much larger than any of the intraband couplings. These arguments are 
supported by the SCDFT result: Using the isotropic approximation we obtain a critical temperature of 8.1 K, against 
T c = 9.4 K obtained within the fully anisotropic approach. Hence, unlike in MgB 2 , we have only a 15% reduction on 
T c in comparison with the anisotropic case. 



V. H UNDER PRESSURE 



The possibility of a superconducting dense molecular phase of hydrogen represents a long-standing problem, recently 
investigated by means of first-principles methods within the SCDFT[llJ. The low temperature and high pressure (> 
400 GPa) phase of molecular hydrogen is predicted to be a base-centered orthorhombic metallic molecular solid 
(known as Cmca phase) [43| with two molecules per unit cell located on different layers. The electronic band structure 
arises from the bonding and anti-bonding combination of the H 2 molecular orbitals. At high pressure, the band 
overlap between the valence and conduction bands produces a complex Fermi surface with disconnected sheets of 
different orbital nature fill]. The molecular nature of the bands provides a strong electron-phonon coupling. The 
superconducting properties of this system are extremely interesting: The calculated T c results to be around 240 K at 
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FIG. 4: Fermi surface of molecular hydrogen at 414 GPa. The different sheets used to construct the interaction matrices are 
labeled as (1): the disk region; (2): the "prism-like" regions; the remaining regions (not labelled). 



~ 450 GPa, with three different gaps on the three different sheets of the Fermi surface. The presence of three gaps 
represents one of the main peculiarities of molecular hydrogen, and its origin is mainly due to the (band) anisotropy of 
the electron-phonon coupling [ll[. However, as we will see, the anisotropy of the Coulomb interaction is fundamental 
to achieve a proper prediction of the gaps at T=0 K. 

In order to simplify the analysis, we will use a multigap BCS model (with parameters calculated ab-initio) to discuss 
the effect of multiband anisotropy, with particular attention to the Coulomb interactions. Thus, we reduce the fully 
anisotropic SCDFT gap equation to a BCS multiband equation, including phonon renormalization effects [which, in 
the SCDFT gap Eq. (fT2"|) are included in the term Z n ^\. In order to do this, we need the band-resolved interaction 
matrices, calculated averaging the k-resolved e-ph and e-e matrix elements over (physically) different sheets of the 
Fermi surface. 

In FigfJJwe show the calculated Fermi surface at 414 GPa. The three gaps separate into three main regions: the disk 
at the Gamma point (labeled as 1), the "prism-like" regions (2) and the remaining. In order to construct a 3-band 
model of the system, we introduce a 3x 3 A- matrix of the partial e-ph coupling \j where the indices i,j = 1, 2, 3 span 
over the above mentioned regions. We also introduce the 3x3 matrix for the V ph , as defined in Sec. IIV CI In analogy 
with the e-ph interaction, we introduce also the 3x3 Coulomb interaction matrix V^ 1 and the matrix defined by 
the relation: /x^ = NjV^. V$ represents the band resolved Coulomb matrix defined in analogy with the definition 
of Xij, but containing the Coulomb matrix elements. In the following, we report the values of these matrices: 



V ph [Ry] : 
2.16 3.00 1.59 
3.00 0.45 0.91 
1.59 0.91 0.78 



V el [Ry] : 
0.71 0.08 0.11 
0.08 0.38 0.20 
0.11 0.20 0.20 







A : 













14 


0.37 


1.45 




0.046 


0.010 


0.100 





20 


0.06 


0.82 




0.006 


0.046 


0.180 
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0.11 


0.71 




0.007 


0.025 


0.182 



We note that the intraband scattering (V? 1 ) is dominating the Coulomb interaction and that it is particularly strong 
for the states of region (1). These states are bonding combinations of anti-bonding molecular states, are responsible 
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for the inter-layer bonding and are very localized in the inter-layer region, with rather poor screening properties. An 
opposite behavior characterizes the \i matrix, which represents the effective e-e interaction in the superconducting 
phase: in fact, despite the large value of V{f we find a very small value of /in, due to the low density of states in region 
(1). Interband terms, in particular those involving the bands less coupled with phonons, dominate thanks to their 
high density of states. Although inspection of the n matrix suggests a minor role of the e-e interaction anisotropy 
in the T c value (we find the same T c in the isotropic and anisotropic case) , we demonstrate that inclusion of band 
anisotropy in the repulsive term has a non trivial effect on the T=0 K gap. The phonon renormalized BCS equation 
(RBCS) is obtained from the SCDFT fully anisotropic equation averaging the k-anisotropy (still retaining the band 
anisotropy) and performing the T — > 0, £, £' — > limits. The results of the solution of the RBCS equation are reported 
in Table I. 





A1/A3 


A2/A3 


Full-RBCS 
RBCS (e-ph) 


1.61 
1.37 


1.13 
1.09 



TABLE I: Results of the RBCS equation including both the e-ph (Xij) and e-e ((J,ij) matrices (Full-RBCS), and without the 
fiij matrix [RBCS (e-ph)]. 

As we can see from Table I the presence of three gaps is predicted even within the RBCS model. Comparing the 
results of the Full-RBCS and RBCS (e-ph), the role of both interaction on the prediction of T=0 gaps is evident. 
Even if the main source of band anisotropy is given by the e-ph interaction, inclusion of the proper band-resolved 
Coulomb interaction increases the anisotropy of about 18%. 

This demonstrates that the band anisotropy of the repulsive e-e interaction, essentially related to the different nature 
of the bands at the Fermi level, increases the multigap superconducting properties in dense molecular hydrogen. 

VI. SUMMARY 

The SCDFT approach allows a fully anisotropic description of the superconducting phase and the inclusion of 
Coulomb repulsion effects on an ab-initio basis. After a brief review of the SCDFT method, in this paper we 
performed a detailed analysis of the Coulomb and electron-phonon matrix elements in intercalated graphite CaC6 
and hydrogen under high pressure. In CaC@ we studied different approximations for the Coulomb interaction. We 
find that, due to the presence of strongly localized states, the use of free electron-like Coulomb matrix elements gives 
a strong underestimation of the Coulomb repulsion and a corresponding overestimation of the critical temperature 
T c . While local field effects are important to describe the screening properties in the graphene layers, inclusion of an 
isotropic Coulomb interaction - averaged on isoenergy surfaces - does not strongly affect the value of T c . Concerning 
the electron-phonon interaction, instead, the isotropic approximation yields a T c reduction of 15%. 

Finally, we calculate the matrix elements of both interactions for H under high pressure. By making use of a 
multiband BCS model using SCDFT-calculated parameters, we demonstrate that the usual isotropic approximation 
of the repulsive e-e interaction is completely unjustified in the case of superconducting molecular hydrogen, and that 
the Coulomb anisotropy it is fundamental for the prediction of the three superconducting gaps at T=0 K. 
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